function [delta_1 delta_2] = func_delta_m_l(q, m, CONSTS)

    k0 = CONSTS.k0;
    a = CONSTS.a;
    eta = CONSTS.eta;
    eta_a = CONSTS.eta_a;
    eps_a = CONSTS.eps_a;
    
    p = func_p(q, CONSTS); 
    [H_m1 H_m2] = func_H_m_l(q, m, CONSTS);
    [J_m1 J_m2 J_m1_1 J_m2_1] = func_J_mk(q, m, CONSTS);
    [n1 n2] = func_nk_from_p(p, CONSTS);
    
    Q = k0*a*q;
    
    delta_1 = -(n2.*((eta/eta_a).*J_m1.*J_m2 - (J_m1+(eta/eta_a)*J_m2).*H_m1) - ...  
              n1.*((eta/eta_a).*J_m1_1.*J_m2_1 - (J_m2_1+(eta/eta_a)*J_m1_1).*H_m1) + ...
              (n2-n1).*(H_m1.^2-((p.^2)./eta_a).*((m^2)./((Q).^4))) + ...
              p.*(eta/eta_a).*(m./((Q).^2)).*(J_m1+J_m1_1-J_m2-J_m2_1)); 
          
    delta_2 = (n2.*((eta/eps_a).*J_m1.*J_m2 - (J_m1+(eta/eps_a)*J_m2).*H_m2) - ...
              n1.*((eta/eps_a).*J_m1_1.*J_m2_1 - (J_m2_1+(eta/eps_a)*J_m1_1).*H_m2) + ...
              (n2-n1).*(H_m2.^2-((p.^2)./eps_a).*((m.^2)./((Q).^4))) + ...
              p.*(eta/eps_a).*(m./((Q).^2)).*(J_m1+J_m1_1-J_m2-J_m2_1));
          
end
